#include <catch2/catch_test_macros.hpp>
#include "InterGridOperator.h"
#include "MultiArrays.h"

TEST_CASE("basic multiarrays funtions", "[MULTIARRAYS]")
{
    MultiArrays<Real, 1> tmp(8, 2);
    auto grid0 = tmp.Level(0);
    for (int i = 0; i <= 8; ++i)
    {
        grid0[i] = i;
    }
    auto grid1 = tmp.Level(1);
    for (int i = 0; i <= 4; ++i)
    {
        grid1[i] = 2 * i;
    }
    // std::cout << tmp;
}

TEST_CASE("interpolation and restriction on 1d", "[OPERATOR]")
{
    MultiArrays<Real, 1> tmp(16, 3);
    GridOp::Linear_Interp<Real, 1> linear_interpolate;
    GridOp::Quadra_Interp<Real, 1> Quadra_interpolate;
    GridOp::FullWeight<Real, 1> FW_Restric;
    GridOp::Injection<Real, 1> Inj_Restric;
    auto grid1 = tmp.Level(1);
    for (int i = 0; i <= 8; ++i)
    {
        grid1[i] = 2 * i;
    }
    linear_interpolate(tmp, 1, true);
    for (int i = 0; i <= 16; ++i)
    {
        REQUIRE(tmp.Level(0)[i] == i);
    }
    Quadra_interpolate(tmp, 1, true);
    for (int i = 0; i <= 16; ++i)
    {
        REQUIRE(tmp.Level(0)[i] == i);
    }
    FW_Restric(tmp, 1);
    for (int i = 0; i <= 4; ++i)
    {
        REQUIRE(tmp.Level(2)[i] == 4 * i);
    }
    Inj_Restric(tmp, 1);
    for (int i = 0; i <= 4; ++i)
    {
        REQUIRE(tmp.Level(2)[i] == 4 * i);
    }
    // std::cout << tmp;

    // TO DO: quadratic interp
}

TEST_CASE("interpolation and restriction on 2d", "[OPERATOR]")
{
    // double h = 1 / 4;
    int n = 16;
    MultiArrays<Real, 2> tmp(n, 3);
    GridOp::Linear_Interp<Real, 2> linear_interpolate;
    GridOp::Quadra_Interp<Real, 2> quadra_interpolate;
    GridOp::FullWeight<Real, 2> FW_Restric;
    GridOp::Injection<Real, 2> Inj_Restric;
    auto grid1 = tmp.Level(1);
    for (int i = 0; i <= n / 2; ++i)
    {
        for (int j = 0; j <= n / 2; ++j)
        {
            grid1[i + (n / 2 + 1) * j] = 4 * i * j;
        }
    }
    linear_interpolate(tmp, 1, true);
    for (int i = 0; i <= n; ++i)
    {
        for (int j = 0; j <= n; ++j)
            REQUIRE(tmp.Level(0)[i + j * (n + 1)] == i * j);
    }
    quadra_interpolate(tmp, 1, true);
    for (int i = 0; i <= n; ++i)
    {
        for (int j = 0; j <= n; ++j)
            REQUIRE(tmp.Level(0)[i + j * (n + 1)] == i * j);
    }
    FW_Restric(tmp, 1);
    for (int i = 0; i <= n / 4; ++i)
    {
        for (int j = 0; j <= n / 4; ++j)
            REQUIRE(tmp.Level(2)[i + (n / 4 + 1) * j] == 16 * i * j);
    }
    Inj_Restric(tmp, 1);
    for (int i = 0; i <= n / 4; ++i)
    {
        for (int j = 0; j <= n / 4; ++j)
            REQUIRE(tmp.Level(2)[i + (n / 4 + 1) * j] == 16 * i * j);
    }
    // std::cout << tmp;

    // TO DO: quadratic interp
}